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We present results of a quenched QCD simulation with overlap fermions on a lattice of volume V = 16 3 x 32 
at (3 = 6.0, which corresponds to a lattice cutoff of ~ 2 GeV and an extension of ~ 1.4 fm. From the two- 
point correlation functions of bilinear operators we extract the pseudoscalar meson masses and the corresponding 
decay constants. From the GMOR relation we determine the chiral condensate and, by using the if-meson 
mass as experimental input, we compute the sum of the strange and average up-down quark masses (m s + rh). 
The needed logarithmic divergent renormalization constant Zs is computed with the RI/MOM non-perturbative 
renormalization technique. Since the overlap preserves chiral symmetry at finite cutoff and volume, no divergent 
quark mass and chiral condensate additive renormalizations are required and the results are O(a) improved. 



1. Introduction 

In the last few years a major breakthrough 
in the lattice regularization of Fermi fields was 
achieved through the closely related domain wall 
|IJ and overlap formulations. Neuberger found 
a lattice Dirac operator D || which avoids dou- 
bling and, most notably, satisfies the Ginsparg- 
Wilson (GW) relation M. Thus the correspond- 
ing action in the massless limit preserves a lattice 
form of chiral symmetry at finite lattice spacing 
and volume ||[| . As a result the use of the over- 
lap action entails many theoretical advantages 
(see ref. || and references therein for recent re- 
views); in particular it forbids mixing among op- 
erators of different chirality and therefore it can 
be very helpful (if not crucial) for computing weak 
amplitudes. No power-divergent subtractions are 
needed to calculate the matrix elements relevant 
for the AI =1/2 rule in K — > 7T7t decays 0. 

The theoretical advantages of overlap fermions 
come at the cost of an increased computational 
burden in numerical simulations. The main ques- 
tion we wish to answer in the work presented here 
is whether the overlap formalism can be effec- 
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tively used for large scale QCD calculations, with 
known algorithms and the current generation of 
computers. The calculation of light quark masses 
is an ideal test since it uses many of the ingredi- 
ents needed for a generic phcnomenological com- 
putation. 

We have performed a fully non-perturbative 
calculation of (m s + rh) in the quenched approx- 
imation following the procedure proposed in ||. 
We also report results for the chiral condensate 
computed using the Gell-Mann-Oakes-Renner 
(GMOR) relation and for the pseudoscalar me- 
son masses and decay constants. 

2. Quark Masses and Chiral Condensate 

The fermionic action in the overlap regulariza- 
tion reads 



S = ip 
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where M is the matrix of bare masses 
(mi, m2,...) in flavor space. The Neuberger- 
Dirac operator is defined as pj 
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where D\y is the Wilson-Dirac operator (defi- 
nitions and conventions used here are fully de- 
scribed in ref. |5)). The renormalized quark mass 
is defined as 



m{jj) = lim Z m (afi)m(a) 



(4) 



where Z m {a\i) is a logarithmically divergent 
renormalization constant and m(a) is the bare 
mass parameter which appears in the action of 
eq. ([j]) and in the corresponding axial and vector 
Ward identities. The non-singlet "local" bilinear 
quark operators we are interested in are defined 



as 



Or{x)=Mx)Y 
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where Or = {V^, A^, S, P} correspond to T = 
{in-, 7m75i 1) 75}- They are subject to multiplica- 
tive renormalization only, i.e. the corresponding 
renormalized operators are 



Ot(x,[i) — lim Zr(a[i)Or(x, a) 



a->0 



(G) 



where Zr{fi) are the appropriate renormalization 
constants. Since V^, and S, P belong to the 
same chiral multiplets Zy = Z4 and Zg = Zp. 
Flavor symmetry imposes Zg = X/Z m . 
The bare chiral condensate is defined as 

X (o) = urn -^(0)[(1 - ^-m\(0)) (7) 

where m in this case is a common mass given to 
the light quarks. For non-zero quark mass 



i(#)[(l-^](0))=x(a) 



(8) 



since chiral symmetry forces the coefficient of the 
linear divergence to be zero. The condensate sat- 
isfies the integrated non-singlet chiral Ward iden- 
tity 



^-$(o)[(i-^M(o)> 

-mV(P(x)P(0)> 
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Therefore by writing the correlation function 
(P(x)P(0)) as a time-ordered product and by in- 



serting a complete set of states in standard fash- 
ion we can also write 
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where Mp is the mass of the pseudoscalar state 
IP). If we use 



2m\(0\P\w)\ = fpMp 



(11) 



where fp is the corresponding pseudoscalar decay 
constant, we arrive to the familiar GMOR rela- 
tion 



X(a) 
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To preserve eq. (|^), the renormalized chiral con- 
densate is defined as 



iV f a^O 

3. Non-perturbative Renormalization 
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The bilinear renormalization constants Zr(afj.) 
have been be computed at one loop in pertur- 
bation theory jlCjJl]] . Nevertheless, to avoid large 
uncertainties due to higher order terms, we prefer 
to compute them non-perturbatively. 

The exact chiral symmetry implies a conserved 
axial current jll],[l2| which enters the correspond- 
ing Ward identities. Therefore the scheme and 
scale independent renormalization constant Za 
can be computed from the axial Ward identity 

Z A {VMx)P{Q)) = 2m(P(x)P(0))+O(a 2 )(14) 

where is the symmetric lattice derivative. 
To compute the logarithmic divergent renor- 
malization constant Zs(fJ-) we implemented the 
RI/MOM non-perturbative renormalization tech- 
nique proposed in ref. |l3f| . The amputated off- 
shell Green's functions are defined as 



A r (p) =S- 1 ( P )G r (p)S- 1 (p) 



(15) 



where Gr{p) are the improved quark correlation 
functions and S^ 1 ^) is the improved (external) 
quark propagator in the Fourier space. The pro- 
jected amputated Green's functions Tr(p) are de- 
fined as 



rr(p) = jTjTr [PtMp)} 



(16) 
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where the trace is over spin and color indices and 
Pr are the Dirac matrices which renders the tree- 
level values of IY(p) equal to unity (see || for 
more details). Z<p(a/z) can be determined, up to 
0(a 2 ), by imposing the renormalization condition 



Zg (afi) = lim Za 



T A (p, m) 



m—*0 



T s (p,m) 



(17) 



To compare the running of Z^ l (afi) with the 
evolution predicted by the renormalization group 
equations, it is useful to define the renormaliza- 
tion group invariant (RGI) renormalization con- 
stant [fill 
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where Cg (//) at the next-to-next-to- leading- 
order (NNLO) is given in [fl5|| . Up to higher order 
terms in continuum perturbation theory Z^ GI (a) 
is independent of the renormalization scheme, of 
the external states and is gauge invariant. It is 
interesting to note the the non-degenerate scalar 
and pseudoscalar Green's functions satisfy the 
WI @ 

(mi - m 2 )T s (p,m 1 ,m 2 ) = (19) 
miTp (p, mi, mi) - m 2 T P (p, m 2 ,m 2 ) 



Equation ( |19[ ) gives an exact relation among 
bare correlation functions of "local" scalar and 
pseudoscalar operators at non-zero quark masses, 
without any reference to the conserved axial or 
vector currents. 

4. Numerical Results 

The numerical results we present are based on 
a set of 54 quenched gauge configurations pro- 
duced with the standard Wilson gauge action, 
with V — 16 3 x 32 and (3 = 6.0, which we 
retrieved from the "Gauge Connection" reposi- 
tory |L6]]. For each configuration we fixed the 
Landau gauge by requiring a quality factor of 
9 < 10~ 6 and calculated the quark propaga- 
tors for 5 different bare quark masses = 
{0.040,0.055,0.070,0.085,0.0100}, as described 
in appendix ^| We computed then the two-point 



correlation functions 
G SS (t) = ^2(S(x,t)S(0,0)) (20) 

X 

Gpp(t) = Y,{P(x,t)P(0,0)) (21) 

X 

GvAp(t) = ^2(V A (x,t)P(0,0)) (22) 

X 

We estimated the errors by a jackknifc procedure, 
blocking the data in groups of three configura- 
tions, and we checked that blocking in groups of 
different size did not produce relevant changes in 
the error estimates. 

4.1. Bare Quark Masses 

In the quenched approximation the contribu- 
tions from chiral zero modes is not suppressed by 
the fermionic determinant. In particular Gpp(t) 
receives unsuppressed contributions proportional 
to 1/m 2 and 1/m, which should vanish in the in- 
finite volume limit but can be quite sizeable at fi- 
nite volume (Tt]] . These quenching artifacts cancel 
in the difference of pseudoscalar and scalar meson 
propagators |Q 



G P ^s(t) = G PP (t)~Gss(t) 



(23) 



The drawback is, of course, that the plateau in 
the effective mass and, correspondingly, the range 
that can be used for the cosh fit become shorter 
(See H for more details). We fitted Gpp(t) and 
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0.089(2) 


0.085 
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0.348(6) 


0.085(2) 


0.070 
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0.315(7) 


0.081(2) 


0.055 


0.0030(5) 


0.280(9) 


0.076(2) 


0.040 


0.0026(5) 


0.239(11) 


0.071(2) 



Table 1 

Meson masses and pseudoscalar matrix elements 
for all the bare quark masses considered in the 
simulations, as obtained from Gp-s(t). 



Gp-s(t) to a single particle propagator in the 
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time intervals t\—t>} — 12—16 and 10— 16 respec- 
tively. The two correlation functions give values 
for the pseudoscalar masses compatible within the 
statistical errors. One does not notice any sign 
of the singular contributions from zero modes in 
the masses obtained from the pseudoscalar cor- 
relation function. These are expected to show 
up at some point, but one would probably need 
much higher statistical accuracy and lower values 
of m to bring them into evidence. The results for 
the matrix elements, parameterised by the factors 
Zs-p and Zpp for Gs~p(t) and Gpp(t), respec- 
tively show more significant differences. Within 
the large statistical errors, Z$-p and Zpp are 
still compatible. However the central values are 
quite different and the fact that a curvature (see 
H for more details) shows up only in the results 
for Zpp points to the fact that what we are see- 
ing is the effect of the unsuppressed zero modes, 
and not of chiral logarithms which would affect 
both sets of results (and would most likely be- 
come noticeable at much smaller values of am). 
On account of the above, we report in table |l| 
the results obtained from Gs—p{t) which we will 
use to derive our further results. It must also be 
said that most of the observables will be calcu- 
lated directly atm~ m s /2 (see below), and for 
these the difference between Zpp and Zs-p is 
irrelevant within statistical errors. We illustrate 
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Figure 1. (aMp) 2 vs. am as obtained from 
Gp-s(t). The dashed lines represent the results 
of a linear fit. 



in fig. |T] the values for (aMp) 2 , obtained from 
Gp-s(t), as a function of the bare quark mass 
am. A linear behaviour 



(aM P ) 2 = A Mp + B Mp (am) 
fits very well the data with 



(24) 



-0.0005(68) B Mp = 1.43(7) (25) 



The vanishing of the intercept within statisti- 
cal errors signals the absence of additive quark 
mass renormalization. In order to fix the lattice 
spacing, we used the method of "lattice physical 
planes" in which the ratio /k/Mk is fixed to 
its experimental value. This avoids recourse to a 
chiral extrapolation for observables except when 
it is really needed. For more details we refer the 
reader to M and only quote here the result 



aM K = 0.216(8) af K = 0.0698(26) 



(26) 



To fix the lattice spacing we compared the value 
of afx with its experimental value and we ob- 
tained aj^ = 2.29(9) We also computed the chi- 
ral condensate by using the GMOR relation. The 
quantity 



a 3 Xm = - (am) 



Z P - 



[aMpf 



(27) 



exhibits a very good linear behaviour in the bare 
quark mass and a linear fit leads to 



a 3 x = -0.00117(27) 



(28) 



The chiral condensate can also be computed di- 
rectly from its definition in eq. (^). Using again 
a computational strategy similar to that used in 
eq. ( p3| ) in order to take care of the infrared di- 
vergent contributions (in m) from unsuppressed 
zero modes, we obtain 



a 3 X (a) = -0.00117(42) 



(29) 



which is in remarkable agreement with eq. (p8|). 

4.2. Renormalization of the Axial Current 

We implemented eq. (|l4|) numerically by com- 
puting the ratio 
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Figure 2. Left: Gsjap/Gpp vs. i for all simulated masses. Right: (ap) as a function of the bare quark 
mass. The dashed line represents the result of a quadratic fit. 



The results obtained for each mass are reported 
in the first plot of fig. The flatness of the 
data highlights the improvement achieved with 
the overlap regularization. We fitted R p {t) to a 
constant p(m) in the range 5 < t < 27 and then 
performed a quadratic fit (second plot in fig. 0) 



p{m) = A p + B p m + C p m 2 
obtaining 



(31) 



1.286(3) (32) 



A p = -0.00002(7) Bp 
C p = 0.277(12) 

Note that the intercept is compatible with zero: 
this should not come as a surprise since A p (x) has 
the correct behaviour under global non-singlet 
chiral transformations. 

The axial current renormalization constant Za 
is given by Za = J- • We also fitted p{m) linearly 
in the quark mass and we take the difference of 
the quadratic and linear results as a rough esti- 
mate of the systematic error. Our final result is 
Z A = 1.55(4) [§. This value is larger than the 
one obtained in refs. JTtJfjJ using standard lattice 
perturbation theory at one-loop. 

4.3. Renormalization of the Scalar Density 

To implement the RI/MOM non-perturbative 
renormalization technique we computed the pro- 
jected, amputated Green's functions of the quark 
bilinears defined in eq. (|l~6|). In the first plot of 
fig. |^ we show the ratio 



Rps = 



(mi - m 2 )T s (p,m 1 ,m 2 ) 
miFp (p, mi, mi) - m 2 T P (p,m 2 ,m 2 



■(33) 



for a given combination of masses. As expected 
from (19) its value is always compatible with one. 
The ratio of the vector and axial Green's func- 
tions produces analogous results which we do not 
report here for lack of space. On the other hand 
the ratio Tp/Ts with degenerate masses is ex- 
pected to be sensitive to non-perturbative Gold- 
stone pole contributions p^ , p^ . A detailed anal- 
ysis will be presented in a forthcoming paper. 

To calculate Zg in the RI/MOM scheme, we 
computed the ratio ||] 

T A {p,m) 



RAs{m) = Z A - , > 
Ts(p,m) 

which we extrapolated to the chiral limit as 



R-Asim) = Z_ 



RI ' Am 2 



Bm 



(34) 



(35) 



to take into account the effect of unsuppressed 
zero modes. The results we thus obtained, shown 
in the second plot of fig. ||, lead to 

Zf\2GeV) = 1.24(5) (36) 

where the error is mainly systematics and is due 
to the uncertainty in the value of the renor- 
malization scale and to the chiral extrapolation 
(the latter estimated by using different forms of 
the extrapolating function) . Using known results 
from NNLO continuum perturbation theory p5|, 
we convert eq. ( |36|) to the MS scheme, obtaining 



Zjf b (2GeV) = 1.41(6) 



(37) 



This result is in good agreement with the one 
in pOj which was determined by fixing the 
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Figure 3. Rps (left) and Zs(afi) (right) as a function of the renormalization scale. 



renormalized quark mass value. Invoking again 
NNLO continuum perturbation theory, we calcu- 
late Zg GI from eq. (|Ts| ). We use Nf — and 
K QCD = 0.238(19) The result, displayed in 
the second plot of fig. || shows the remarkable 
flatness of Z^ GI for > 1.5 GeV which con- 
firms the high level of improvement reached by 
the overlap regularization. Averaging over the 
region (1.5 — 3) GeV, we find 



Z 



RGI 



0.82(3) 



(38) 



5. Physical Results 



In this section we will use our lattice results to 
infer the renormalized values of the sum of the 
strange and average up-down quark masses and 
of the chiral condensate. From eq. (Eg) we obtain 



m s (a) + m(a) = 149(9) MeV 



(39) 



where the error is statistical only. Since our vol- 
ume is fairly large, we expect our main sources of 
systematic errors to come from discretization ef- 
fects and from the quenched approximation. For 
a rough estimate of the systematic error due to 
quenching approximation we can use the results 
in p^] , from which one sees that, within the 
quenched approximation, the use of alternative 
observables to calculate the lattice spacing can 
produce differences of approximately 10%. Had 
we used tq |^3| to fix our lattice spacing we would 
have obtained a number ~ 7% higher than the 
one in eq. (B9|). Combining this fact with the 
results in |24|], we infer that discretization uncer- 
tainties are below 10%. If we had used the ex- 



perimental NNLO results from a s (Mz) — 0.118 
in the matching of the renormalization constants 
in eq. (|3^), we would have obtained a value of 
the scalar renormalization constant ~ 10% higher 
than the one given above. The difference can 
be taken as an indication of the systematic er- 
ror introduced by the quenched approximation. 
In order to be conservative, we will take 15% as 
the estimate of our overall systematic error in the 
renormalized quark masses due to quenching and 
discretization effects. A more precise estimate of 
the systematic errors will need much more ex- 
tensive simulations, which would go beyond the 
capability of our current computer resources and 
the exploratory nature of the our work. Combin- 
ing the results in eqs. ( |36| ) and (|3^) we obtain 

(m s + to) ri (2 GeV) = 120 ± 7 ± 21 MeV (40) 

which represent one of the main results of this 
work. By using the quark mass ratio m s /m = 
24.4±1.5 from chiral perturbation theory [£5[ the 
above translates to 



MS 



(2 GeV) = 102 ± 6 ± 18 MeV 



(41) 



This result agrees very well with the current lat- 
tice world average 

Insofar as the value of the condensate is con- 
cerned, if we used the standard two-step ap- 
proach, i.e. first measure the dimensionless con- 
densate, see eqs. (^) or (^9|), and then multiply it 
by the cubic power of the lattice spacing, the re- 
sult would be affected by a very large systematic 
error due to the uncertainty in the determination 
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of the lattice spacing in quenched simulations. In- 
stead, we will use an alternative method [E7| ■ We 



write the GMOR relation (12) for the renormal- 
ized condensate as follows 



X(a) 



(42) 



where Bm p is defined in eq. (E5J) and f x = 
0.1282 GeV is the "experimental" value, in phys- 
ical units, of the pseudoscalar decay constant ex- 
trapolated to the chiral limit. While computing 
the condensate from the above formula relies on 
an additional element of experimental informa- 
tion, it has several advantages. The most im- 
portant is that, by expressing the condensate in 
terms of / x , we are left with only one power of 
the UV cutoff a" 1 . With this method we obtain 

(^) m (2 GeV) = (43) 
- 0.0167 ± 0.0010 ± 0.0029 GeV 3 

where the estimate of the systematic error has 
been made using the same criteria we used to 
compute the error for the quark masses. This 
is our best value for the chiral condensate. It is 
interesting to note that, if we had used the stan- 
dard technique, starting from eq. (p8f), we would 
have obtained a result with a central value very 
close to the value in eq. (^3|), but with a much 
higher error. Finally, from eq. ( fill ) and NNLO 
matching, we get 

(tpipf^p GeV) = - (267 ± 5 ± 15 MeV) 3 (44) 

This result is in very good agreement with the re- 
sult obtained by the authors of refs. |28|]29] 1 , while 
it is smaller than the result in [^0[ , even if still 
compatible within errors. Our result is also com- 
patible within errors with the number obtained 
few years ago in p7| with Wilson-type fermions. 
We expect, though, the systematics due to the 
discretization effects to be smaller (0(a 2 )) in the 
result reported in eq. (Q) than the error (0(a)) 
which affects the determination in [^7f . 

6. Conclusions 

At this conference we presented results for the 
pseudoscalar meson masses and decay constants 
together with the first fully non-perturbative 



computation of (m s +m) with overlap fermions in 
the quenched approximation. We also computed 
the chiral condensate (^ip) from the GMOR 
relation and directly. To avoid uncertainties 
due to lattice perturbation theory, we com- 
puted the multiplicative renormalization con- 
stant Zs(/mi) non-perturbatively in the RI/MOM 
scheme. While the systematics errors due to 
quenching are common to previous calculations, 
the other systematic errors (mostly discretization 
effects) arc different than in other lattice regular- 
izations and likely to be smaller, because of chiral 
symmetry. Our results have indeed produced a 
remarkable verification of "good chiral behavior" 
of the overlap fermions both in the axial Ward 
identity and for the pseudoscalar masses. 

The calculation of light quark masses uses 
many of the ingredients needed for a lattice cal- 
culation of weak matrix elements, although the 
latter is computationally more demanding. From 
this point of view, the very good agreement be- 
tween our results for the quark masses and the 
current lattice world average bodes well for the 
use of the overlap formalism also in matrix el- 
ement calculations. Our investigation has been 
mostly of exploratory nature. One would need 
to extend it to larger volumes and better statis- 
tics. Nevertheless, we believe that it demon- 
strated that the overlap formalism can be used ef- 
fectively, with known algorithms and the present 
generation of computers, for large scale QCD cal- 
culations, at least in the quenched approxima- 
tion. Thus we would conclude that it represents 
a very promising non-perturbative regularization 
for solving long standing problems, such as the 
proof of the AI —1/2 rule and the calculation of 
e'/e. 

A. Implementation of the Overlap Opera- 
tor 

The complicated form of Neuberger operator 
renders its numerical implementation more in- 
volved and expensive than those of the most com- 
mon rcgularizations. The crucial point is to im- 
plement the sign function 



e(H) = sign(tf) 



(45) 
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over the whole spectrum of eigenvalues of H = 
J5X. Exact diagonalization is feasible j3l[] for 
two dimensional models, but becomes prohibitive 
for larger systems such as QCD. Many algorithms 
have been proposed in the literature f3^-|34|l . We 
opted for the optimal rational function approxi- 
mation suggested in refs. f32| , |34|| . Starting from 
the observation Q that the sign function can be 
written as 



s(x) 
s N (x) 



lim e 7v (x) 

N^oo 



X Cq + 



N 



Ck 



qk 



(46) 



a good approximation can be obtained at finite 
N by fixing the coefficients c k and q k with the 
Remez algorithm p4] , p5[ . 

The numerical procedure we have followed to 
obtain a quark propagator from a fixed source 
vector 1 77), i.e. solving the equation 



(l__L mi )£>. 

2p 



is the following: 



\Xi 



\v) 



(47) 



• We have extracted the 15 smallest eigen- 
values of H 2 and the corresponding eigen- 
vectors minimizing the Ritz functional ]36| , 
projected them out and treated the action 
of e(H) on this subspace exactly. The mini- 
mum precision required for each eigenvalue 
A, 2 is 



fritz 



\H 2 \\ t ) - X 2 ^ 2 



(A, I A, 



< 10 



-6 



(48) 



where |Ai) is the corresponding eigenvector. 
We checked on all the configurations that 
the largest eigenvalue extracted is always 

\Xraax\ > 0.15. 

• We have rescaled the reduced matrix H = 
H/r by a factor r — 3.7 and we have ap- 
proximated e(H) as in eq. ( [46| ) with N = 14 
and Ck,qk fixed with the Remez algorithm 
p4lB5[ which guarantees a precision 



\e(H) - e N (H)\ < 3.6 x 10" 



(49) 



in the range H e [0.040,1. 
invert all the equations 



H 2 + q k )\ipk) = \s) 



N 



In order to 



14 (50) 



in one stroke we use an (inner) multi con- 
jugate gradient (MCG) solver @. Note, 
that during the inner MCG it is necessary 
to store only N + 2 large vectors since we 
are only interested in 



JV 



(51) 



k=l 



The convergence is governed by the inver- 
sion corresponding to the smallest q k for 
which we require a residual 



|(ff a + «H)hk 



\s)\ 2 



< 10~ 6 (52) 



Notice that having extracted the lowest 
eigenvectors of H 2 improves the approxima- 
tion used for e{H) and reduces the condi- 
tion number of H 2 , which results in a speed 
up of the inversion in eq. (|5C 



• From eq. (g 
write 



and using the GWR we can 



\Xi) 



2p 



)Dt 



(p 2 - m 2 /4) (D + £>t) + pm ; 



■h)(53) 



and therefore we are left to invert 



D + D 



pmf 



p 2 — mf /4 



m = \v) (54) 



D + D< is an Hermitian matrix and there- 
fore a MCG can be applied also in this 
case. Moreover [D + D + ,^ 5 ] = implies 
that we can always use sources \rj) and so- 
lutions \xi) restricted to one chiral sector, 
saving one application of s(H) per itera- 
tion J34j] . The convergence is governed by 
the inversion corresponding to the smallest 
quark mass, for which we require a residual 
of r nut < 10~ 5 
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We did not attempt to use any preconditioning 
for the inner and the outer MCG. At the i-th step 
of the outer MCG, the difference \di) of two con- 
secutive source vectors |s<) and |s»_i) on which 
en(H) has to be applied gets smaller and smaller 
when Tout decreases. Therefore at the i + 1 step 
we have applied en{H) directly on \ ck) with an in- 
creased r in (every 40 outer MCG step we required 
always the full precision to avoid rounding accu- 
mulation). After convergence, we always checked 
that the true residual associated with the eq. (47) 
for each mass is < 10~ 5 . In our simulations the 
average numbers of inner and outer iterations is 
0(150) and 0(250) respectively. 
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